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We study the mixed state in an extreme type-II lattice dj.2_y2-wave superconductor in the exper- 
imentally most relevant regime of intermediate magnetic fields Hd <^ H <^ Hc2- We analyze the 
low energy spectrum of the problem dominated by nodal Dirac-like quasiparticles with momenta 
near ~ (±fc_D, ifcu) and find that the spectrum exhibits characteristic oscillatory behavior with 
respect to the product of ko and magnetic length I. The Simon-Lee scaling, predicted in this regime, 
is satisfied only on average, with the magnitude of the oscillatory part of the spectrum displaying 
the same dependence as its monotonous "envelope" part. In general, the spectrum obeys a 

scaling law E„ii — — p-£„(kZ, t/A, knl), where f is a dimensionless universal 27r-periodic function of 

kol- The oscillatory behavior of the spectrum is due to the inter-nodal interference enhanced by the 
singular nature of the low energy eigenfunctions near vortices. Our results constitute an example 
of a finite size scaling of the Dirac-type quantum criticality. We also study a separate problem of 
a single vortex piercing an isolated superconducting grain of size L x L. Here we find that the 
periodicity of the quasiparticle energy oscillations with respect to koL is doubled relative to the 
case where the field is zero and the vortex is absent, both such oscillatory behaviors being present 
at the leading order in . Finally, we review the overall features of the tunneling conductance 
experiments in YBCO and BSCCO, and suggest an interpretation of the peaks at 5 — 20 meV ob- 
served in the tunneling local density of states in these materials. We find that in the case of a pure 
d-wave superconducting order parameter with featureless vortex cores, the zero bias conductance 
peak (ZBCP) appears only on the sites that are the immediate nearest neighbors of vortex locations, 
while all the other sites in the close vicinity of vortices exhibit no such ZBCP, and instead display 
pronounced peaks at sub-gap energies, typically at about a half or less of the coherence peak energy. 
Furthermore, we find that the on-site ZBCP can be strongly suppressed by enhanced local pairing 
near a vortex. 



PACS numbers: 74.25.Jb, 74.25. Qt 
I. INTRODUCTION 



The study of the quasiparticle spectrum in the mixed 
state of d-wave superconductors followed sooiiSi^ af- 
ter the d-wave nature of pairing in the cuprate super- 
conductors became apparent. A distinct character of 
the Tunneling Local Density of States (TLDOS) in s- 
wave and d-wave superconductors near a vortex was pro- 
posed by Wang and MacDonald to serve as a test which 
would allow determination of the pairing symmetry in 
the cuprateai. For d-wave pairing, they found that the 
TLDOS at the vortex location plotted as a function of en- 
ergy, exhibits a prominent peak at zero applied bias volt- 
age, while in an s-wave case, for otherwise similar param- 
eters of the model, the thermally broadened TLDOS has 
a minimum at E — surrounded by two large sub-gap 
peaks. While their calculations apply to densely packed 
vortices and unrealistically high magnetic fields in real 
cuprates, a similar conclusion concerning this "zero-bias 
peak" is obtained also within a single vortex calculation 
of Franz and Tesanovic^. Interestingly, the STM exper- 
iments in BSCCO^ and YBCO^ - now unambiguously 
known to be of a d-wave type - reveal that the "zero- 
bias peak" is completely absent. Instead, the tunneling 
conductance experiences a dip at zero bias, and new rel- 
atively small sub-gap peaks at energies 5 — 20 meV. To 



explain this discrepancy - which at the moment is unre- 
solved - one usually relies on additional order parame- 
ters in the vortex cores. This line of thought advocates 
that due to the suppression of the superconducting order 
parameter within vortex cores, a new competing (local) 
order emerges there, which ultimately is responsible for 
the deviations from Refs. m. Several different order pa- 
rameters were considered in the literature: d + id super- 
conducting order—, antiferromagnetic order—, pseudogap 
state^, circulating current&W'^^, d-density-wavei^. Other 
explanations of the absent ZBCP include the anisotropy 
of the tunneling matrix elementsi^, and, most recently, 
quantum zero-point motion of vorticesi^. 

In addition to the above "high-energy, short-distance" 
features of a vortex core - an interesting problem in its 
own right - the quest for a description of nodal quasi- 
particles within some simple low-energy effective Hamil- 
tonian, which would facilitate theoretical analysis, was 
launched as a separate line of inquiry. An effective de- 
scription is important in order to address more compli- 
cated problems such as the interactions of fluctuating 
vortices with nodal d-wave quasiparticles or the effects 
of disorder in a nodal superconductor. The initial steps 
in this direction were taken by Simon and Lee^^ who pro- 
posed that, after extracting the rapid "kp" oscillations 
of the wavefunctions, the "linearized" effective version of 
the Bogoliubov-deGennes (BdG) Hamiltonian suggests a 
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simple scaling ()15|l for the low energy {E ^ A) sector 
of the quasiparticle spectrum in the mixed state of type 
Il-superconductors, and consequently for various other 
measurable quantities. The scaling function was then 
calculated by Franz and Tesanovic-S, who employed a 
singular gauge transformation (FT transformation) and 
expansion of the wavefunctions in the plane wave ba- 
sis to find that the spectrum at the very low energies 
E ^ hvp/l is essentially the same as the original spec- 
trum of the zero-field problem, but for the renormaliza- 
tion of the slopes of the anisotropic Dirac cones at the 
nodes. The linearized FT Hamiltonian was subsequently 
analyzed both numerically a nd the oretically, using its 
symmetry properties, in Refs. llTllSi 

Further study, however, revealed several new ques- 
tions. Quite separately from its origin, the linearized 
Hamiltonian turned out to be somewhat challenging to 
analyze due to singularities at vortex positions which ren- 
dered it incomplete unless its proper self-adjoint exten- 
sion is constructed by imposing an additional boundary 
condition at each vortesi^Si^. Such boundary conditions, 
which turned out to be necessary in performing the nu- 
merical and symmetry analysis of the linearized Hamil- 
tonian, are discussed at length in a companion paper^^. 

Furthermore, the relation of the linearized description 
to the tight-binding model also turned out to be some- 
what more complex than initially anticipated: the Simon- 
Lee scaling of the quasiparticles energy eigenstates ac- 
cording to ((T^ demands that if the spectrum is gap- 
less on the linearized level - as found in Ref. - the 
gaps of the full non-linearized problem must decrease as 
l/l'^ as a function of I or faster in the limit of small 
magnetic fields. However, exact diagonalization of the 
non-linearized problem at zero chemical potential /i = 0, 
showedi^ that the gap in the spectrum oscillates between 
and 0(hvp /I), depending on the commensuration of the 
magnetic length to the atomic lattice spacing. Perhaps 
the most telling manifestation of the intricate relation be- 
tween the linearized continuum and the tight-binding lat- 
tice models of the mixed state is the exact resulfe^S for the 
spectrum of the latter when = and 1/5 = 2 (mod 4): 
in this case the number of the zero energy states is dou- 
bled compared to the zero magnetic field result. Clearly, 
such doubling is difficult to account for if one uses the 
non-perturbed plane wave basis as the departure point 
for a perturbation theory. 

Here, we explore further the non-perturbative effects of 
the tight-binding (TB) model of a d-wave superconduc- 
tor in the presence of a vortex lattice. In section H we 
present the results of a systematic study of the spectrum 
for large magnetic lengths I (low magnetic fields, corre- 
sponding to realistic values in cuprates), up to Z = 120(5, 
where 6 is the lattice spacing, and for general /i. We 
start by focusing on the low-energy properties of the spec- 
trum and analyze the validity Simon-Lee scaling for this 
model. The dispersion is shown to obey the scaling on 
average^ and in general to experience rapid oscillations of 
the energy levels as a function of both the magnetic field 



and /i. These anomalous oscillations, which can be uni- 
fied within a new generalized form of Simon-Lee scaling, 
are described by an additional dependence of the energy 
levels on the commensuration of the inter-nodal distance 
and magnetic length I. The inadequacy of the Simon-Lee 
scaling in its conventional form is shown to be the result 
of the singular nature of the BdG eigenfunctions com- 
bined with the inter-nodal interference, as conjectured in 
Refs. Liaad The linearized effective Hamiltonian is ar- 
gued to still accurately represent the low energy sector 
of the theory, but the necessary condition is stricter than 
anticipated earlier and demands also /cf^ S> 1 rather 
than only kpl 3> 1. 

In section HI we describe the high-energy, short- 
distance features of the spectrum. We find that although 
the TLDOS is indeed peaked at the four sites of the tight- 
binding lattice surrounding the vortex, in agreement with 
previous work, the four immediate neighbors are rather 
an exception than a rule: all sites in the proximity of 
the vortex - except the nearest neighbors - exhibit no 
zero-bias peak, and furthermore have additional peaks 
at sub-gap energies. In the concluding section, we dis- 
cuss how the on-site peak can be suppressed and argue 
that the 5 — 7 meV peaks observed in STM experiments 
could in fact be due to regular c?-wave vortices, but with 
a particular profile for the amplitude of d-wave gap func- 
tion on those few bonds constituting the cores. 

The anomalous enhancement of the inter-nodal inter- 
ference by singular potential due to vortices has also a 
prominent effect on a related single vortex problem, i.e., 
an isolated superconducting grain of size L x L in a com- 
mensurate magnetic field H = $o^~^j where the elemen- 
tary flux $0 equals hc/{2e), pierced by a single vortex. 
This is discussed in detail in section IV, where we show 
that although the spectrum has oscillations of the low 
energy levels of magnitude proportional even in zero 
magnetic field due to finite-size effects, in the presence of 
the vortex the periodicity of these oscillations is doubled. 

II. TIGHT-BINDING LATTICE MODEL 

A. BdG equations 

We start from the Bogoliubov-deGennes(BdG) Hamil- 
tonian Htb of the modeiiSiSS, which is defined by its ac- 
tion on a two-component Bogoliubov-Nambu wavefunc- 
tion 7/>r = {urT^r)^ as 




In the simplest case, the hopping and pairing fields de- 
scribed by trr' and Arr' are nonzero only on the nearest 
neighbor bonds, and in the presence of magnetic field B, 
the hopping irr' ~ t*,^ is modified by Peierls factors 

trr' = —texp{—iArr') = —texp i / A • dl ) , (2) 
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where A is the vector potential corresponding to mag- 
netic field B. The pairing field/gap function Arr' should 
in principle be determined from a self-consistent proce- 
dure stemming from the same microscopic Hamiltonian 
that in the mean- field approximation yielded . For ex- 
ample, in the simplest model that results in the d-wave 
order within the mean-field approximation for a wide 
range of parameters - the extended Hubbard model with 
the nearest neighbors density-density attraction grij-nri - 
the self-consistency condition reads 



En 



Arr' = g ^(WrV + ""r'^r ) tanh ^ 



(3) 



where T is the temperature and {umVn) are the eigen- 
states of the BdG Hamiltonian of energy En- While in- 
corporating the self-consistency condition is not an im- 
possible task, the results to a certain extent will depend 
on the microscopic model, from which the condition was 
derived. 

In the context of the superconductivity in cuprates, 
however, such dependence is very weak: the amplitude 
of the order parameter |Arr'| recovers rapidly to its uni- 
form state value A, while the phase is subject to the 
condition of overall winding by 2-k along any lattice path 
enclosing a vortex. These two conditions suggest a useful 
simple Ansatz for the starting point of the iterative self- 
consistency procedure: Arr' = Aryrr' exp(i0rr')j where 
the bond phase 6'rr' is given in the Appendix b and the 
d-wave nature of the bond field Arr' enters through fac- 
tors ryr,r+5 = 1 {'rir,r+s = -1) if ^ = ±5x {5 = ±6y) . 
One then proceeds to diagonalize Htb from Q, recom- 
putes Arr' using the self-consistency condition ||2J), and 
repeats the procedure until the convergence is achieved. 
In practice, the starting Ansatz is a very good approxi- 
mation to the final solution in a sense that both have the 
same phase defects, the same symmetries, and the ratio 
of the fully self-consistent solution Arr' to the Ansatz 
A77rr' exp(i0rr') IS merely a periodic smooth function 
close to unity at all bonds of the lattice, except possibly in 
a close proximity of vortices. Consequently, we first con- 
centrate on the (non-selfconsistent) BdG Hamiltonian 
with Arr' = Aryrr' exp(i0rr')7 which allows for an easier 
systematization of the results due to the reduced num- 
ber of parameters and note that all possible microscopic 
Hamiltonians leading to c?-wave pairing are now encoded 
by a single bond variable A. At the same time, this ap- 
proach permits one to avoid a time-consuming search for 
the solution of the fully self-consistent problem. In the 
end, we will briefly discuss the effects of varying Arr' 
near vortices. 



vortices, Hamiltonian Htb does not explicitly possess 
the symmetries of the physical vortex lattice. In gen- 
eral, one has to accompany the translations by a vortex 
lattice vector with an additional gauge transformation 
restoring the Hamiltonian to its original form. Rather 
than working with representations of the resulting mag- 
netic translations group, we perform a unitary transfor- 



mation of a special form U = diag(e* 



) such that 



the transformed Hamiltonian H = U^^HtbU is explic- 
itly periodic. 

It is easy to see that, regardless of the transformation 
used to bring the Hamiltonian Htb to a periodic form, 
the minimal unit cell must contain at least two vortices: 
after the transformation the diagonal hopping term con- 
tains modified Peierls factors exp(iylrr') = exp(i^rr' + 
iar' — ictr). Suppose these factors are periodic, then con- 
sider a product n{rr') e*'*"'"'' over the oriented bonds along 
a closed path formed by the perimeter of the unit cell tra- 



versed counterclockwise. Since e 



and fac- 



tors exp(zArr') are to be periodic, such a product should 
be equal to 1. 



(rr'> 



{rr'> 



le 
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Thus the flux of magnetic field through the unit cell must 
be an integer of l-Kh/zj e — hc/e^ and therefore must con- 
tain at least two hc/2e superconducting vortices. 
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FIG. 1: Left: magnetic unit cell containing two vortices la- 
belled as A and B with magnetic length / = 85. Right: The 
symmetry of dispersion within the Brillouin zone (BZ), 
which follow from the symmetry operations of the Hamilto- 
nian H, is shown. 16 equivalent points are displayed as solid 
dots; and it is sufficient to study only 1/16-th portion of the 
BZ drawn as a dashed triangle. 

An explicit form of such a transformation can be real- 
ized by considering a simple family of the so-called sym- 
metric transformations22i^: 



H 



■0r' - fJ-OSIpr 



.(4) 

where the summation indices r' denote the nearest neigh- 
bors of r and is the Pauli matrix. Although we will 
be interested in a periodic inversion-symmetric lattice of 



U, 
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diag(e^"',e-'"0, 



with suitably chosen a-^. In the continuum version of 
the theory this transformation forces branch-cuts and 
non-locality on the eigenfunctions of -ff^^. Although 
the tight-binding lattice version of it does not cause un- 
due complicationSiS^, here we utilize another common 
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choice, the FT transformationi&, whose continuum ana- 
logue does not require branch-cuts: 



U = diag(e''^- , e 



(5) 



where site variables (j)^ and 4>b are explicitly given in the 
Appendix A of Rcf. 21. After defining periodioi^iSl bond 
variables V^, , V^, , and Arr' according to 



(6) 



Ar' =err' -0^'-</'^?, (7) 

the resulting periodic Hamiltonian is given by 




with r' denoting the nearest neighbors of r. 

B. Linearization of the tight-binding Hamiltonian 

One can now attempt to describe the low energy por- 
tion of the quasiparticle spectrum using the linearized 
approximation, leading to Simon-Lee scaling for various 
properties. The derivation of the linearized version of 
Htb has been performed in Refs. 19 20 by gradient ex- 
pansion resulting in a continuum Hamiltonian with dis- 
persion, which is then linearized as usuali&. One might 
expect that the coefficients of the final model obtained 
in this way reproduce the spectrum of the full TB model 
for the values of chemical potential fi not too close to the 
half-filling - where the dispersion is not quadratic ~ and 
also for /i not too close to the bottom of the tight-binding 
band, otherwise the linearization procedure itself is not 
justified. Below, we describe an alternative derivation 
which, while proceeding along similar lines as the stan- 
dard procedure, does lift the first restriction and allows 
one to consider values of jjL at and near half-filling. At 
the end of this section, we will revisit this derivation and 
show that it should be corrected to accommodate the 
curvature and the so-called "interference" effects. 

The linearization procedure is based on the assumption 
that wavefunctions can be represented as 



E 



(9) 



where j labels the Dirac-like nodes of a d-wave gap func- 
tion located in momentum space 

and -0'^-'^(r) is slowly varying function on the scale of 
kp^. Variable kjj introduced here for brevity of nota- 
tion simply equals where kp is the magnitude of 
the Fermi momentum in a nodal direction. After sub- 
stituting this form into the BdG eigenvalue problem for 
Hamiltonian _ff (jS)) , a typical term has the form of a sum 
over 5 = ±(5i;, ±(Sy: 
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S 



r+S 



(10) 



where Fourier transform of /r is assumed to be composed 
of wavevectors close to the four nodal momenta at the 
Fermi surface and V • w = 0. Note that the straightfor- 
ward gradient expansion is valid only qualitatively since 
kpS is not necessarily much smaller than one. The lin- 
earization, however, can be performed directly, without 
the preliminary "continuization" step, by first separating 
the rapidly oscillating part: 



fr 



where is a function that changes slowly on the scale 
of few lattice spacings. To obtain the effective linearized 
description, we now replace the lattice function F^ by a 
slowly changing function F(r) defined in continuum, such 
that it coincides with F^ when r correspond to the lattice 
sites. Thus, F{r) can be thought of as an interpolating 
function for a discrete set F^. The detailed form of the 
interpolation turns out unimportant for the leading order 
results derived below, as long as the characteristic scale 
on which F^ varies is much larger than the lattice spacing 
S. Then, in the expression 



S = 



s 



we use 



and expand slowly varying function F{r) into Taylor se- 
ries while retaining factors exp(ik^'' • S). The result for 
node k^^-* = {kajkn) is 



2V2Ssm{kDS) 



Px +Py 



F 



V2 V2 

-6^ COs{kDS){p + ^^)^F +...], (11) 



where p denotes the usual continuum momentum opera- 
tor — iV, and . . . denote terms containing higher powers 
of (5V and Sw. 

The expressions for the off-diagonal terms of ^ differ 
only by the presence of the factor rjs: 



and after similar algebra we find 



S' 



2V2Ssm{kDS) 



Py - Px , Wy -Wx 



F 



V2 V2 
- cos{kDS)((p^ + w^)^ - {py + Wyf)F + ... , (12) 

Using the expressions for S and 5', in the leading order 
we obtain 

Hun = VF y=-^(T3+VA ;^ (Ji+Vp (13) 
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V2 
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where the effective velocities vp and va are given by the 
zero-field expressions given in the next subsection (|18|l . 
V is the superfluid velocity, and Hi — pi + ai is the gen- 
eralized momentum with a = (v^ — vb)/2 describing 
the vector potential due to an array of Aharonov-Bohm 
TT-fluxes located at vortices of subset A and similarly — tt- 
fluxes at vortices of subset B. Since the only length in 
this effective low-energy Hamiltonian is the magnetic 
length I, it immediately follows 

HHnir,l,VF,VA) = ^^Hii,i{r/1,VF/VA), (14) 

and, consequently, the spectrum of this Hamiltonian 
must satisfy the Simon-Lee scaling relation for the low 
energy eigenstates: 

E^{k)^^E:{kl,VF/vA) ■ (15) 

C. Zero field spectrum 

Let us briefly summarize the properties of the spec- 
trum in the absence of a magnetic field, which is obtained 
from ||SJ) by setting the field-induced factors cxp{iArr') 
and exp(zV^l^') to zero. The excitation spectrum in this 
case is 

ek = ±^J^l + Al, -tt/S < k^^y < tt/S 

where = —2t{cos -t- coskyS^ fi and the d-wave 
pairing gap function is Ak = 2A(cosfc2;^ — cosfc^^). The 
four nodal points (i = 1, . . .,4) of the spectrum are lo- 
cated at k^^ = (±fc£),±fc£)), where 

fc_D = i arccos (-|^) ■ (16) 

All four nodes merge at fi — ±4i, while at /i the 
inter- nodal separation is the largest. The dispersion in 
the vicinity of each node can be approximated as ek — 
^VpSk"^ + ^^(Jfc^, where Sk±{Sk\^) is the displacement 

of the momentum from a nodal point in the direction 
perpendicular (parallel) to the Fermi surface, and the 
effective velocities are 




D. Spectrum in a magnetic field: /j. — (known 
results) 

The structure of the spectrum in a finite magnetic field 
is a great deal more complex. In all cases, the spectrum 



at energies much larger than A evolves to a set of sharp - 
conventional Schrodinger, as opposed to Dirac - Landau 
levels, while at low energies is characterized by strongly 
dispersive energy bands. For concreteness, we will con- 
sider a square lattice of vortices oriented as shown in 
Fig. n with a unit cell of minimal area, which contains 
two vortices, labelled A and B. In addition, we will only 
study the values of magnetic field that correspond to the 
symmetric placement of vortices within plaquettes - this 
requirement restricts the values of magnetic length I to 
even integers in the units of lattice spacing. 

Let us first recall the results for a fully particle-hole 
symmetric systemi^iSSiSS, corresponding to /i = 0. Due 
to a special symmetry (ipr (— l)*^^+y-'/'^?/'r) of this case, 
the spectrum is doubly degenerate for all momenta k. 
As shown in Ref. 123 . when the center of inversion for 
the vortex lattice coincides with a site of the atomic lat- 
tice - a situation realized for magnetic lengths 1/5=2 
(mod 4)- the spectrum contains eight Dirac nodes {six- 
teen zero energy states): two degenerate nodes at each 
of the four momenta k = (±|^,±|^). This is quite un- 
usual non-perturbative result, since it suggests that for 
arbitrarily small fields giving rise to a vortex lattice, the 
number of zero modes is doubled compared to the four 
nodes of the zero-field problem, provided the magnetic 
length has the correct commensuration with the tight- 
binding lattice spacing. 

In the opposite case 1/5 = Q (mod 4), the symmetry 
does not demand the zero modes, and, as was found nu- 
merically in Ref. ITfll the field-induced gap in this case 
scales as as a function of magnetic length within the 
regime Hd <^ H <^ Hc2- This result is also surprising, 
since for the lowest energies, one expects to recover the 
Simon-Lee scaling of the linearized problem (|15|l . 

The existence of the nodes for commensuration 1/5=2 
(mod 4) might suggest that the spectrum of the lin- 
earized problem given by E* in (|15|l is gapless - this is 
also a result obtained earlier^^--'^^ from the analysis of the 
linearized Hamiltonian. This conclusion in turn would 
have required that in the expansion of the overall field- 
induced gap 

1 1 
Am = aij + ai— + . . . 

the leading l/l term is absent, and only small gaps of 
order l/l'^ should generally appear from the terms that 
were left out in the process of linearization. 

The large gaps whose size scales as l/l for weak mag- 
netic fields in a particle-hole symmetric situations at 
commensu ration factors 1/6 = (mod 4), as was found 
m Refs. Il9l2fl therefore come entirely unexpected. 
These large gaps were interpreted as as the effect of the 
inter-nodal interference, and such effects for the delib- 
erately distorted lattice were indeed found to be sup- 
pressed. It was argued that in realistic situations a weak 
disorder in vortex positions or the one due to impuri- 
ties will suppresses these interference effects. Yet, the 
following questions remain to be answered: First, is this 
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situation specific only to a particle-hole symmetric sys- 
tem (/z = 0)? Second, how could the scahng relation l(T3|) 
be explicitly violated, even for an ideal periodic lattice? 
Finally, how - if at all - the Simon-Lee scaling can be re- 
covered in the tight-binding problem without introducing 
the disorder explicitly? 



E. The spectrum for general /i 

We start by answering the first of these question and 
consider the spectrum for non-particle-hole symmetric 
systems (/i 7^ 0). We show below that the rather in- 
volved behavior displayed by the spectrum as a function 
of /i and I in fact follows a simple universal pattern when 
expressed in terms of suitably rescaled variables. For 
the square lattice of vortices considered here, the anal- 
ysis is further simplified due to the sixteen-fold symme- 
try of the dispersion En]^ within the Brillouin zone illus- 
trated in Fig. ^ Furthermore, using a transformation 
■0r (— l)(^+^)/*'i/)r and the symmetry of the spectrum 
at each k as a function of energy (see Appendix) , it is easy 
to show that the spectra at values of chemical potential 
+fi and —fi are identical, with or without magnetic field, 
and therefore in what follows we assume fi > 0. 
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FIG. 2: Upper panel: The overall gap Am, induced by vortex 
lattice as a function of chemical potential fi. The parameters 
are /S. = t, 1/5 = 18, 22, 30 (solid symbols) and / = 365 (open 
symbols). Center panel: the gaps are renormalized by l/vp 
according to the Simon- Lee prescription: Am — Aml/fiVF- In 
the limit of small magnetic field the result should have been 
independent of in direct contradiction with explicit numeri- 
cal evaluation shown here. Lower panel: instead of the chem- 
ical potential /i, the horizontal axis represents knl/i'iiT). In 
the limit of low magnetic fields (i S> 5) all curves representing 
dependence of Am on Iko collapse onto a 27r-periodic func- 
tion. For fixed I, deviations from this universal scaling are the 
largest for fi close to the bottom of the tight-binding band, 
where the Fermi surface is small and the validity condition 
for linearization {Ikp S> 2n) is violated. 
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FIG. 3: Upper panel: The dependence of the lowest eight pos- 
itive energy bands E„k at fixed momentum koZ = (7r/2,7r/2) 
on chemical potential fi. The parameters are A = t, l/S = 
30, 36. Note that the two values of magnetic length I corre- 
spond to the gapped {l/S = 4n) and gapless {l/S = in + 2) 
families at half-filling (/i = 0). Center panel: the energy lev- 
els at ko are rescaled by l/vp. Lower panel: Instead of the 
chemical potential /i, the horizontal axis represents lkD/{2n). 
In the limit of low magnetic fields {I 2> S) all curves repre- 
senting dependence of Am onkol collapse onto a 27r-periodic 
function. 

The field- induced gap for vp — va is plotted as a func- 
tion of the chemical potential /i in the upper panel of 
Fig. |5| for several values of magnetic length I. At all 
magnetic fields the dependence on fj. displays a charac- 
teristic oscillatory behavior. For the family of magnetic 
lengths I = (4n + 2)S, the spectrum is gapless at /i = 0, 
in accordance with the previous resultsiS, and for a finite 
fraction of one cycle of oscillations in /i the field-induced 
gaps are extremely small, their size quite possibly set by 
in the scaling limit: we were not able to definitely 
establish the scaling behavior due to the smallness of the 
gaps in this regime. Then the gaps increase and remain 
large - of order Hvp/l - for about a third of the cycle, 
until eventually again turning to zero. This cycle is then 
repeated over and over. For I ~ AnS the only difference 
is that the cycle is offset by half a period in /i. The over- 
all slow decrease of the average gap size for large values 
of /i follows directly from the Simon-Lee scaling (|15|l as 
Vp decreases with fi (see Eq. (|18|l '). To account for the 
expected Simon-Lee scaling, the central panel of Fig. |21 
shows the rescaled gap 



as a function of fi. If Simon-Lee scaling in its original 
form were exact, one would expect to see no dependence 
fi. Instead, for any given field value, the rescaled gap 
itself exhibits oscillatory behavior. 

Still, comparing the upper panel of Fig. |21 one must 
conclude that is a step forward compared to A,„; on 
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average the curves representing different magnetic field 
and difi^erent values of fi look almost identical. Although 
the detailed analytic theory of the "interference effects" 
remains a challenge for the future, the essence of such 
interference is vividly illustrated by replotting the family 
of Am as function of Ikn /{"^Tr). So rescaled, all curves 
collapse into a single universal periodic function shown 
in the third panel of Fig. |21 

We find that the above oscillatory behavior is not spe- 
cific to the field-induced gap function; the dependence 
of the entire spectrum i5„(Zk) is characterized by similar 
behavior. As an example, Fig. |21 displays the eight low- 
est energy levels at k = (0, 0) for I — SOS and I = 36(5, 
representing two families of magnetic fields. The central 
panel shows the energy levels rescaled by vp/l, while the 
bottom panel shows that the remaining oscillations of the 
rescaled spectrum fall onto a universal periodic curve if 
plotted as functions oi kjjl rather than ^. 




13 14 15 16 17 18 19 20 21 



lk^/(2Tl) 

FIG. 4: Upper panel: The dependence of the lowest eight 
positive energy bands _E„k at fixed momentum koZ = (0, 0) 
on chemical potential fj,. The parameters are A = Q.2t, I — 
505. Lower panel: Instead of the chemical potential fi, the 
horizontal axis represents Zfco/(27r). Only the range < /i < 
3.5t is plotted in the lower panel. 

The pattern just described is not restricted to the 
isotropic case a a = vf/va — 1 and holds for all a a 
we checked (2, 5, 10, 20, 50). For large anisotropics, the 
deviations from the scaling behavior become more pro- 
nounced at smaller values of / due to violation of the 
hiOc <C A condition, which translates to l/S 3> \/aA - 
this condition is necessary if we are to treat magnetic field 
as a "small" perturbation of a zero-field d-wave supercon- 
ductor. Fig. ^illustrates the spectrum at fc = for the 
anisotropic case q;a — 5. In all cases, for sufficiently weak 
magnetic fields, the spectrum still exhibits 27r-periodic 
oscillations in k]jl of amplitude Chvp/l, where the nu- 
merical prefactor C is only a function of t/A. We stress 
that the oscillatory part of the spectrum is of the same or- 
der of magnitude as its smooth, "envelope" dependence. 



and is not smaller than this average Simon-Lee envelope 
part in any sense. These results can be summarized in 
the following scaling form: 

Snk = ■^f„(k?,t/A,A:i30, (19) 

where £ is a universal dimensionless function, which is 
27r-periodic with respect to its last argument. This scal- 
ing, which combines the oscillations with respect to both 
magnetic field and the chemical potential, holds for all 
values of chemical potential, except when /i is very close 
to the bottom of the band - in this case the scaling was 
studied by Vafek and Tesanovic^^. 

The origin of this oscillatory behavior can be traced 
back to the linearization procedure. Recall that the 
matrix elements of the FT-transformed, but yet non- 
linearized Hamiltonian 7i (|SJ) are evaluated in the plane 
wave basis. By inspection, in the limit Ik^ 3> 1, the 
leading term of this infinite matrix Tikk' appears to have 
a block-diagonal form, with each of the four blocks cor- 
responding to separate nodes. Within each block, only 
the leading order approximation in is kept, while 
the block-offdiagonal (inter-nodal) matrix elements as 
well as the higher-order corrections to the block-diagonal 
(intra-nodal) matrix elements, which superficially scale 
at worst as l^^ , are neglected. Therefore, by construc- 
tion, the matrix elements of each block precisely coin- 
cide with the matrix elements of the linearized Hamilto- 
nian. What we found, howeveriSSiSi, is that the situation 
is not this simple: this description is necessarily incom- 
plete due to the singular character of the linearized prob- 
lem wavefunctions near vortices. Consequently, the lin- 
earized Hamiltonian requires self-adjoint extensions, ob- 
tained by imposing additional boundary conditions near 
vortices, which themselves are ultimately determined by 
those "higher-order terms" which were dropped in the 
oversimplified analysis. 

A straightforward way to appreciate the significance 
of these subdominant terms is actually to recompute 
the matrix elements of the full Hamiltonian - not in 
the plane-wave basis as one naturally would within a 
pedestrian perturbation theory - but with respect to 
the exact eigenf unctions of the linearized problem^^ with 
fixed boundary conditions at vortex locations {9}. The 
wavefunctions of the linearized Hamiltonian diverge near 
vortices^-'- as r"^/^, and therefore the eigenfunctions of 
the linearized Hamiltonian at each node can be written 
as 

where i labels the nodes 1, 1, 2, and 2, k is the Bloch mo- 
mentum, n is the band index, and dimensionless continu- 
ous functions x and / are such that the expression in the 
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brackets is Bloch-periodic. Note that if the wavefunctions 
contained only the regular part, the terms retained in 
the linearized Hamiltonian would be of the order hvp/l, 
while the non-linear terms such as etc, 

would contain an additional factor of ln(//^), 
and could have been safely omitted. The presence of 
the divergent part of the wavefunctions described by the 
first part of (|20l) . however, changes the situation. Let 
us evaluate again the structure of the matrix elements 

(^nk''l'^l*n'k') bctwccn the states at momenta k and 
k' differing by reciprocal lattice vector G. The intra- 
nodal matrix elements with j — j' to the leading or- 
der are just the eigenstates of the linearized Hamilto- 
nian hvpE^{\<d)5nn' /l- The corrections due to non-linear 
terms, however, quite peculiarly also exhibit the same 
scaling as a function of magnetic length I as we will ar- 
gue now. Consider a typical non-linear term since 
(mv)^ increases near vortices as down to distances 
of the core size '--^ ^, to the leading order the matrix ele- 
ment 



cx / {riy^/^ — {rl)-^/'^{rdr) 



oc 



Vp 1 



(21) 



has the l/l scaling. Higher-order curvature terms con- 
taining higher order derivative operators and potentials a 
or V can be estimated similarly; corrections to the matrix 
elements due each successive term \ai^ Y^dl^''^ in general 
are of the order of VF{kF£,)~''-'~^^ /I- Therefore, the lin- 
earization procedure in the presence of magnetic field is 
justified when the condition {kpC) 3> 1 is satisfied, and 
the role of the small parameter is played by both {kpC)^^ 
and {kpl)~^ <^ 1, and not only the latter, as is commonly 
assumed. 

Moreover, the "interference" terms relating different 
nodes (j 7^ j') have a similar form and scale as 



0) 



(/) 



Vp 1 



le 



2e 



(22) 



where G = (k^'^ + k) - (k^"^ 



k') and coefficients 



C1.2, determined by the wavefunctions Xk "^ ^-nd 
(q = A,B), depend on (ri, k;n',k') but not on k^^* or 



One therefore generally anticipates that - with 



kp^ kept fixed, as in our model, and other parameters 
(such as I or fi) freely varied - the spectrum will undergo 
a complicated evolution, even at the leading order in 
(c.f. Ref. .IQi) . due to the inter-nodal contribution en- 
hanced by the singular character of wavefunctions near 
vortices. Note that in the tight-binding lattice model 
with the nearest neighbor hopping terms only, we are 
precisely in this situation: in our simplified model, where 
no self-consistency condition is employed, kp^ is a fixed 



number of order 1 since the role of the cut-off ^ is played 
by the lattice spacing S, and kp is bounded by {tt/2)S~^. 
Even when the self-consistency condition is employed, as 
long as the uniform system is described by nearest neigh- 
bors only, the ratio vp/v^ automatically fixes t/A; for 
a fixed anisotropy of d-wave nodes, kp^ is bounded by a 
number of the order of Q!a since vp/v^ t/ A kp^ and 
£, S, and therefore the simple Simon-Lee scaling limit 
will be difficult to reach in the strict sense. Of course, 
one may introduce the next-nearest neighbors and fur- 
ther hopping terms in order to optimize parameters and 
maximize kp£^ while retaining the fixed value of ckA- In 
this case, the amplitude of the oscillations will still scale 
as l"^, albeit with a suitably reduced prefactor. 

The above "interference" pattern of the spectrum for 
a moderately large (kp^) is expected to have a periodic 
structure, depending on the commensuration of the nodal 
wavevectors and a magnetic length. Consider a change 
in chemical potential /i or other parameters that result 
in a displacement of nodal points at the Fermi surface. 

If the difference (k^'' — k^ ■*) • (R^ — R^) changes by 
a multiple of 27r, then the amplitudes of the matrix el- 
ements in (|22ll between (n,k) and (n',k'), which deter- 
mine the leading order perturbative corrections to the 
energy spectrum, are the same to the leading order in 
apart from the prefactor {kp£,)~^. Thus, in addition 
to overall the Simon-Lee %p/r' scaling, the spectrum has 
periodic oscillations determined by the commensuration 

(1) d') 

of the inter-nodal momentum Icp — kf, and the differ- 
ence Ra — Rs- More precisely, the spectra for two sets 
of parameters will be similar whenever the nodal points 
{zLkp),ztkp)) and (ztkp)' ,zLkp)') satisfy the condition 

kol-ko'l' ^2TTn, (23) 

where n is an integer. This is equivalent to Eq. H19|) 
surmised from the numerical solution. 

A remarkable feature of the oscillations seen in Figs. 
121 and 121 is that the frequency of oscillations in fi grows 
rapidly with magnetic field. Incidentally, this suggests a 
way of incorporating the effect of weak disorder, which is 
expected to suppress the oscillations, in a relatively sim- 
ple manner: weak disorder in a full calculation is equiv- 
alent to weakly modulated /x(r). On the other hand, 
since the spectrum is a rapidly oscillating function of fi, 
only the quantities averaged over one period of oscilla- 
tions are of interest. For a typical value of magnetic field 
(~ 1 Tesla) the period of oscillations can be estimated 
to be of order 10 meV. Thus, if a random impurity po- 
tential /i(r) is of comparable magnitude or larger, only 
the averages of measurable physical quantities over a pe- 
riod of oscillation are observable. In Fig. [3 the density 
of states (DOS) averaged over the first and the second 
periods close to half-filling are shown for a variety of 
magnetic fields. At low energies, DOS is linear in en- 
ergy and should be associated with the nodal structure 
of the spectrum on average. 

We end this section by alerting the reader to the fact 
that, although we have performed detailed numerical cal- 
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FIG. 5: Full DOS for q:a ~ 1 averaged over one cycle of os- 
cillation in fj, for several values of magnetic field characterized 
by magnetic length The density of states was computed for 
all values of p on a mesh of size Sfi = 0.02t for I = 225, 305 
and Sn = O.Olt for I — 365. The averaged DOS is shown for 
the first period of oscillations p G (0, 4t sin(27r5//)) and for 
the second period ^ £ (4t sin(27r5/Z), 4f sin(47r5/Z)). The inset 
shows the enlarged low-energy part of the figure. 



culations for the nearest-neighbor hopping model only, 
where the strict reduction to the simple linearized de- 
scription is hampered by the enhanced effects of the inter- 
nodal interference and curvature terms, we expect that in 
a more elaborate (and more realistic) model, with longer 
range hoppings, where the condition kp^ 3> 1 is bet- 
ter satisfied, the linearized effective theory does indeed 
provide a quantitatively faithful description of the low 
energy sector of the theory. In that case, the singular na- 
ture of the potential due to vortices still requires special 
care, but such care can be administered by constructing 
suitable self-adjoint extensions of the linearized Hamilto- 
nian, as belabored in Ref. |^ where a detailed discussion 
of the Dirac-Bogoliubov-deGennes quasiparticles in sin- 
gular vortex potentials is presented. 



III. TLDOS NEAR VORTICES AND THE 
MISSING ZERO-BIAS COHERENCE PEAK 

So far we were describing the details of the spectrum 
at the lowest energy scale set by hvp/l- Now we turn to 
large-scale properties of the TLDOS g{r, E) - a quantity 
of direct interest in the STM experiments, which can be 
expressed through the eigenstates of the BdG Hamilto- 
nian (wr, Wr) as 

g(r,i?) cx ^(|7.„k(r)|V'(£^- £^„k) 

nk 

+ k„k(r)|V'(£^ + £^nk)), (24) 




E/t E/t 



FIG. 6: Left panel: TLDOS at = in a mixed state at 
4 representative points for I — 505 (lines) and I = 385 (sym- 
bols). Far from the vortices, the TLDOS of a uniform d-wave 
superconductor is recovered. At four corners of the plaque- 
tte containing the vortex the TLDOS exhibits the zero-bias 
peak', while at the next-nearest and the next-next-nearest 
sites, the TLDOS develops peaks at sub-gap energies. Not 
only the position of these peaks, but also the thermally broad- 
ened TLDOS at all energies does not depend on magnetic field 
(length) provided that the temperature is larger than a typi- 
cal width of a band (T — 0.05t here). This large-scale, "high- 
energy" behavior of the thermally broadened TLDOS should 
be contrasted with the stark dependence of the low-energy 
features on I, commensuration effects etc, which we focused 
on in the previous section. This "fine" structure, which corre- 
sponds to true TLDOS is shown on the right, where TLDOS 
is plotted for I = 225. As we described earlier, the latter 
is generically gapped with the gap scaling as l^^ or, as in 
the example shown in the right panel, is linear for special 
commensuration between the magnetic length and the Fermi 
momentum. 



where r denotes the site of the TB lattice and E is the 
bias. While the results at these large energies are less 
universal and depend to a much larger degree on the 
band structure, the spatial profile of the order parameter 
etc., certain qualitative features turn out to be rather 
robust and will be discussed in this section. 

A typical dependence of the TLDOS on bias for 
particle-hole symmetric case fi — is shown in Fig. El 
for a set of representative points of the tight-binding lat- 
tice. First note that the thermally broadened TLDOS 
is essentially field-independent once the temperature ex- 
ceeds the typical width of the field-induced bands, which 
varies from C{a/^)hvF /I at low energies, where the spec- 
trum is strongly dispersive to huc ~ ^irt/P at energies 
larger than Ac (see Refs. IT^It^IT^ . As an example, 
compare the TLDOS for / = 50(5 (Hues) and I = 30(5 
(symbols), which are shown in the left panel of Fig. I^jfor 
temperature T = 0.05i. 

Far from the vortices the TLDOS is similar to the zero 
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magnetic field result, and quite naturally the deviations 
grow progressively stronger as one approaches a vortex. 
As found by Wang and MacDonald^, on four cites sur- 
rounding the plaquette with the vortex the thermally 
broadened TLDOS has a pronounced maximum at _E = 
as a function of the applied bias voltage, which is called 
the zero bias conductance peak (ZBCP). Note, however, 
that ZBCP appears only after thermally broadening the 
TLDOS, which in its original form is either gapped at 
general /i or has a linear dependence on the energy for 
a discrete set of fj, such as = and 1/6 = 2 (mod 4) 
as discussed in the previous section (see the right panel 
of Fig. To resolve this low-energy gapped or linear 
dependence, however, the temperature must be smaller 
than a typical width and separation between the bands. 
We stress that the ZBCP does not correspond to any "lo- 
calized state" : many narrow Bloch bands (see right panel 
of Fig. El merge into a peak after thermal broadening. 

Importantly, the situation is very different on the next 
nearest and next-next nearest neighbor sites around a 
vortex: there the local density of states exhibits peaks at 
energies w ±Ac/2, where Ac denotes the coherence peak 
energy in a uniform system. Note that these peaks share 
several similarities with the sub-gap features observed in 
experiments, namely, the energy of these sub-gap peaks 
is independent of magnetic field (see left panel in Fig. ^ 
and it also increases with Ac. Again, these peaks do not 
correspond to any "localized state(s)", as many narrow 
bands contribute to the peaks after thermal broadening 
of the LDOS. Although in itself this observation does not 
quite suffice to explain the absence of ZBCP in experi- 
ments, it does suggest that the experimentally observed 
TLDOS might be reproduced quantitatively by consid- 
ering a standard d-wave vortex on a lattice with some 
relatively minor modification, rather than invoking the 
appearance of additional order parameter(s) within vor- 
tex core(s). 

A hint of such a minor modification, which could sup- 
presses the ZBCP at the four sites closest to the vortex, 
comes from a recent analysis of the dopant oxygen atoms 
in BSCCO. As was noticed by Nunner et ali^, the nature 
of spatial correlations between the position of the oxygen 
atomsSS and features observed in the TLDOS, indicates 
that the strength of the electron-electron effective pair- 
ing coupling constant, and therefore also the magnitude 
of the d-wave superconducting gap function in BSCCO, 
are both strongly susceptible to local variations. Vari- 
ations of A by a factor of two or more on a scale of a 
few lattice spacings, which are basically never seen in 
conventional superconductors, are routinely observed in 
BSCCO and other cuprates. In a zero-field case it was 
found'^^ that these modulations of the A are likely to 
be caused by dopant atoms; while the detailed micro- 
scopic origin is not clear at this point, it is conceivable 
that dopant atoms cause local distortions of the atomic 
lattice and cause spatial variation of the superexchange 
interaction or other interaction important for supercon- 
ductivity. Since vortices are expected to be pinned by 



(0,0) 

(0,1) 

(1,1) 




FIG. 7: Left: the thermally broadened TLDOS for a d- 
wave superconductor with spatially constant amplitude of the 
gap function, except on four bonds surrounding each vor- 
tex, where the gap function is set to zero. The parameters 
are: A = 0.2t, I = 385, ^ = 0. The temperature is cho- 
sen as T = 0.05t, and the origin (0, 0) denotes the upper- 
right corner of a plaquette containing the vortex. Note that 
the nearest sites to the vortex exhibit the ZBCP, while sites 
(1,0) and (1,1) have sub-gap (sg) peaks at non-zero energies 
O.SAag. Compared to the case where there is no suppression 
of |Arr'| near vortex core, the ZBCP is slightly enhanced. 
Right: Same, but bond variables on four bonds around each 
vortex are increased by a factor of 3. Note that this change 
affects most dramatically the nearest sites (0, 0) to the vortex, 
where the ZBCP is suppressed and two maxima at sub-gap 
energies develop. 



impurities, the profile of the order parameter near a vor- 
tex would consequently differ from that in an ideal model 
considered so far. Furthermore, even if the vortex is not 
pinned by an impurity, it may distort the lattice and 
cause variations of the order parameter near its core that 
are not described in the canonical BdG scheme, where 
the pairing interaction constant is assumed to be spa- 
tially uniform. 

While at this stage the above reasoning in the context 
of cuprates is only a speculation, it is still a useful phe- 
nomenology to examine more closely the effect of such 
modulation of the d-wave gap function near the vortex 
core. The main results are summarized in Fig. [7| First, 
consider a suppression of the gap function all the way 
to zero on four bonds surrounding plaquette with the 
vortex; the result is shown in the left panel of Fig. |7| 
The ZBCP is strongly enhanced, while other features are 
modified only a little. In the opposite case (see right 
panel of Fig. \7^, when the magnitude of d-wave gap func- 
tion is locally enhanced, the zero-bias conductance peak 
is strongly suppressed, the spectral weight is transferred 
to the Ac/2 sub-gap states, and TLDOS acquires a form 
rather similar to that observed in experiments. The sup- 
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pression of the ZBCP is even stronger if the bonds with 
enhanced gap function extend further. Note that in ex- 
periments the sub-gap peaks have an energy lower than 
Ac/2: in BSCCO Pan et al. reported^ the sub-gap(sg) 
features at Agg = 7 meV for the samples with the co- 
herence peak at Ac = 32 meV, while Hoogenboom et 
alm^ cite values Agg — ±14 meV and Ac = ±45 meV. 
In the earliest samples of YBCO, where these low-bias 
features were observed, Maggio-Aprile et al. reported 
Asg = ±5.5 meV for system described by Ac — ±20 — 25 
meV. In all cases, the ratio Agg/A^ ranges from 0.20 
to 0.33. Since we used the simplest tight-binding model 
described by only two parameters {t,fj,), the numerical 
discrepancy between the result Asg/ Ac ~ 0.5 and the 
experimentally observed ratios is not unexpected. 

Finally, we comment on the 4x4 modulations observed 
in the vicinity of vortex cores^Si^. Such modulations are 
likely to be caused by strong fluctuations of the d-wave 
order parameter, which are believed to become enhanced 
near vortices. Explanation of such effects is clearly be- 
yond the scope of the present paper based on a simple 
mean-field formulation. It has, however, been argued by 
several groups^SiSiS^iS^ that even in the absence of mag- 
netic field the enhanced phase fluctuations of the d-wave 
gap function result in a broken translational symmetry 
with modulated local average of the gap function. This 
provides an alternative mechanism of the modulations in 
the absolute value of A. 



IV. SINGLE VORTEX PROBLEM 

In this section, we study a problem of nodal BdG quasi- 
particles in presence of a single vortex piercing a super- 
conducting grain or droplet of size Lx L, where L relates 
to an external magnetic field in such a way that exactly 
one superconducting flux $0 = hc/{2e) fits into the sys- 
tem. A continuum version of a similar problem was con- 
sidered in Ref. fH, where the delocalized character of the 
core quasiparticle states was established. 

The present problem is technically somewhat easier to 
handle than the periodic array of vortices of the previous 
sections; however, there are certain general features com- 
mon to both situations. In particular, the anomalous en- 
hancement of the inter-nodal interference and curvature 
terms by the r~^/^ increase of the wavefunctions near 
vortex location within the linearized framework still in- 
fluence the spectrum, albeit now in a less dramatic fash- 
ion - unlike the translationally-invariant case considered 
in previous sections, quasiparticles within an isolated su- 
perconducting grain have energy levels that to the leading 
order in exhibit oscillations as function of kijL even 
in the absence of magnetic field. When a magnetic field 
is turned on, the singularities of the wavefunctions near 
the vortex result in the halving of the oscillation period. 

The starting point for description of the quasiparticles 
inside such a superconducting grain is still the Hamilto- 
nian except now the BdG wavefunctions {u{r),v{r)) 



are required to vanish outside the grain. Alternatively, 
all bond variables such as Apr' or trr' can be set to zero 
on links along the perimeter of the grain. The remaining 
bond variables Arr' in principle should be determined 
from the self-consistency conditions, however, just as in 
the case of the vortex lattice problem, this approach has 
a drawback of depending on the precise form of the mi- 
croscopic theory and furthermore on the precise nature of 
the boundaries. Following our justification from the pre- 
vious section, and in order to focus on the simplest model 
with the least number of parameters, we describe the re- 
sults for the order parameter with a constant amplitude 
I Apr' I, its phase simply given by the polar angle around 
the origin. We verified explicitly that after implement- 
ing the self-consistent solution of the problem using the 
condition (|3Jl we find only small quantitative deviations 
from the results described in this section, and no change 
in the qualitative conclusions. Although the "constant 
amplitude" , "polar angle" approximation for the order 
parameter is violated near the boundaries of the grain, 
its effect on the physics near the vortex appears insignif- 
icant. 

Before presenting the numerical results, let us start 
with several simple observations. Consider again the low- 
energy effective description of HBdG- The singular gauge 
transformation and subsequent linearization proceed just 
as in the case of the vortex lattice problem with the re- 
sult given by (|13|) . which can be rewritten in the scaling 
form p4l) . A significant difference at the level of the 
linearized description is the presence of external bound- 
ary. Although the rescaled Hamiltonian H'^^ contains 
no information on the system size L and the microscopic 
lengthscale kp, the spectrum still does not exhibit the 
scaling 

En = ^FniaA) , (25) 

where Fn is a universal function of the anisotropy a a, as 
one might initially suspect. Instead, both the boundary 
conditions at the grain's edge and the singular character 
of the wavefunctions near the vortex affect the leading or- 
der result for the spectrum of the Hun, and consequently 
violate the simple scaling relation 12511 . whose more ap- 
propriate form should be 

En = ^Fn{aA,B.C.) . (26) 

The label B.C. here stands for boundary conditions deter- 
mined by dimensionless combination ikjjL)^ that naively 
might have seemed to drop out of the leading part of the 
scaling. 

As an illustration, consider first a simple example of a 
zero-field problem - a lattice d-wave superconductor in an 
empty box with impenetrable walls. The eigenfunctions 
of the non-linearized problem are given by 

V'fe..fc„(r) = Csm{k,x)sm{kyy) W . (27) 
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FIG. 8: Left: Superconducting island of size L x L. The 
tight-binding lattice is shown in black, the red circle is the 
vortex, and the dashed bonds correspond to the boundary 
of the superconducting region with both t^yi and A^r' = 
across the dashed bonds, solid lines. Right: Four nodes of a 
d-wave superconductor in momentum space at (±fci3,±fc_D), 
where ko = arccos(— /i/(4f)), for a general incommen- 
surate case are shown as four circles. The horizontal and 
vertical lines display the grid of commensurate wavevectors 
{■KUx/L, imy/L). 



where C is a normalization constant, and due to the zero 
boundary conditions at the edges, we have ki = imi/L 
with positive integer Ux and Uy. 

The components of the Nambu spinor Uk and t;k can 
be expressed as 

= (1 + ek/Sk)/2 , 

vi = {i- ek/i?k)/2 , 

where = Sk + ^k- ^^o^c that the solution ^ mixes 
four plane waves exp(zk • r) with k — (±kx,iky) in a 
specific combination, and no other combinations are al- 
lowed. Among other things, it suggests that simply tak- 
ing eigenfunctions of the linearized Hamiltonians corre- 
sponding to energy E and combining them in all possible 
combinations in ip^-^^ Q will not work: only special su- 
perpositions, which make the full wavefunction vanish at 
the edges of the system, are allowed. Consider now the 
lowest energy levels: if the node 1 situated at {ku,k]:i), 
coincides with one of the allowed mesh points in mo- 
mentum space {mix/ L,Trny/L), then the lowest energy 
value is simply zero. Otherwise, depending on anisotropy 
au the lowest energy level is reached at one of the four 
points of the mesh closest to the node {k^jki)). More 
precisely, if /cd = nn/ L + Sk with \6k\ < 7r/(2L), then the 
ground state energy is given by the least of ET^n/L^irn/L 
and £^7rn/L,7r(n+sgn(5fc))/L, which to the leading order in 
1/L equal 



E. 



■Kn/ L,-Kn/ L 



A\tsm{kD5){5k)\ 



(28) 



and 



E. 



Therefore, the ground state energy is given by H28|l when 
\Sk\ < (1 -t- AVi^)7r/(4i), and by ^ otherwise. Note 
that the result is an oscillating function of kjj changing 
from zero, whenever kjjL — nn and the nodes coincide 
with mesh points, to 



21 sin 



m{kDS)\^t^ (2\Sk\ - + (^j ' 



(29) 



2n / 

27rmin(i, A)sin(fcD(5)/i = min(i, A) — W 1 - . 

L V loi^ 

As a result, strictly speaking, there is no uniform scaling 
(|25|l of Simon-Lee type even for the zero-field problem, 
no matter how large the magnetic length L is. Instead, 
the scaling is fulfilled only on average. 

Why do the linearization, and the scaling relation l|25|) 
fail? The answer is that in this problem the Fermi mo- 
mentum scale kp has not truly been eliminated from the 
linearized problem. Although the linearized Hamiltonian 
^lin does not contain any dependence onkp, the bound- 
ary conditions that should be satisfied by the eigenfunc- 
tions of the Hamiltonian retain the information on com- 
mensuration between ki) and 1/L. To derive them in 
general, consider again ©. At the left boundary x = 0, 
and the condition reads 

+ e'^'^^y (V'(^) (0, y) + ^® (0, y)) = . (30) 
Since ^pi vary slowly on the scale of 1/kp, we obtain 

V'W(0,y)+V'''H0,y) = , (31) 
V'(^^(0,y)+V''^(0,y) = 0. (32) 

Similarly, the boundary condition at the right edge of the 
square is 

e'^'^y (e^'^^^V^^) (i, v) + e-'^"^i/^^ {L, y)) 

(33) 

and consequently 

g»fcoL^(i)(2,^ y) + e-'''"^i;<^^\L, y) = (34) 

e-^'^^^V^^n^, y) + e'''"^^^^^{L, y) = . (35) 

In addition, similar conditions must be satisfied at y = 
and at y = L. Note that the conditions on the 
eigenfunctions of the linearized problem couple different 
nodes, and moreover they explicitly involve a phase factor 
eicp{2ik]jL). It is expected, therefore, that the spectrum 
of the problem does depend on kiyL modulo tt even in 
the leading order of approximation. 

The vortex problem adds a new layer of complexity 
when the linearization is performed. Although in the full, 
non-linearized problem, the divergence is cut-off at about 
vortex core radius ^, and does not pose complications, 
the effects of the curvature terms such as the inter-nodal 
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FIG. 9: Left panel: The energy of the lowest 20 energy levels as a function of n is shown in black solid circles. The three panes 
of the figure correspond to A = f/4, A — t/2, and A = t. In all cases, the size of the system L equals 41. The red dashed lines 
are the energies of the lowest four eigenstates of the zero-field problem. Center panel: Rescaled version of the left panel. The 
rescaled energy E„L/{hvF) levels are plotted versus {ko — n /{2S))L/2tt. The 7r/(25) offset is artificially introduced to make 
the plots with different L collapse on the same graph near half filling /x = 0. Right: The energy eigenvalues rescaled as in the 
center panel are shown for L — 415 and L = 1215. Only the first cycle of the oscillatory pattern is shown for compactness. In 
both cases, A = t. 



interference are enhanced due to the singular character of 
the vortex potential at the linearized Hamiltonian level, 
as explained in the previous section. A simple estimate 
similar to H22() shows that as a result, the oscillations of 
the energy levels are controlled by condition H23|) . and 
thereby the periodicity of the energy levels in a grain 
vifith a vortex, plotted as a function of k^L, should be 
doubled compared to the oscillations in an empty grain. 

The spectrum of a tight-binding lattice superconductor 
of size L X L found numerically in a presence of a single 
vortex in a magnetic field HL^ = hc/{2e), is shown in 
Fig. El (left panel) . The low-energy quasiparticle spec- 
trum exhibits the following properties: (i) a generic spec- 
trum is gapped, (ii) for A < Ac, where Ac « 0.75i, the 
zero energy states appear at discrete values of the chem- 
ical potential /i, (iii) the spectrum does not follow a sim- 
ple scaling relation H25() displaying instead an oscillatory 
behavior, with the magnitude of the oscillations scaling 
as hvp/L, and the period of oscillations decreasing away 
from /i = 0. The last property is a direct consequence 
of the approximate 27r-periodicity of the spectrum with 
respect to koL: the spectra therefore must be similar 
at /i and /x' related by L{kD{p) — kD{pL')) = 27rn, where 
koifJ-) — arccos(— /i/(4f), and therefore the periodicity 
condition for / 3> 6, when the equivalent values of /i are 
closely spaced, can be written as 

Just as in the case of a vortex lattice, it is useful to redis- 
play the data by extracting the analogue of the overall 
Simon-Lee scaling factor vp/L from the energy levels, by 
plotting EL/vp vs koL (see Fig. center). Clearly, 
after such rescaling, the dependence of the energy levels 
on /i (or ko) is more uniform compared to the raw data 



plotted in Fig. |5| The oscillatory part of the spectrum, 
however, also scales as 1/L (see Fig. O right panel). 
Note that the periodicity of the oscillations in koL is 
27r, twice the periodicity of the zero-field problem, as ex- 
pected. The pattern of Fig. El holds extremely well for 
all ^, except near the extreme values close to /i « — 4t, 
where the kp becomes comparable with 1/(5, and the lin- 
earization procedure is not justified. The commensura- 
tion effects could only be experimentally accessible at the 
temperatures smaller than the amplitude of the oscilla- 
tions in energy (« hvp/L for the lowest energy levels) 
The thermally broadened quantities such as the TLDOS 
will not reflect these oscillations unless the temperature 
T ^ hvp/L. Besides the temperature, impurities, in- 
strumental resolution and other factors could result in 
broadening the energy levels. 

Now we turn to large-energy, short-distance features 
of the quasiparticle spectrum and describe the TLDOS 
calculated at temperatures larger than the separation be- 
tween the energy levels, but still smaller than A. This 
implies sufficiently large L. In practice we used typically 
T = 0.05t and L ranging from 30(5 to 120(5. A represen- 
tative TLDOS is shown in Fig. EH The fcuZ-oscillations 
at the lowest energy scales C{a^)hvp/L are essentially 
absent in such thermally broadened LDOS. While the re- 
sulting TLDOS will still depends on /i (or kp^l)^ the de- 
pendence is not oscillatory and merely reflects the slow 
varying changes of the normal state band structure, yield- 
ing large-scale changes such as the position of the van 
Hove peak. Overall, the spacial and energy distribution 
of TLDOS is quite similar to the vortex lattice case: the 
TLDOS at the center of the vortex has a zero-bias co- 
herence peak, while far from a vortex and the edges the 
TLDOS is similar to the uniform zero- field result, as ex- 
pected. At the 4 nearest and 4 next nearest neighbors 
the TLDOS has pronounced peaks at 1/3 — 1/2 of the 
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FIG. 10: TLDOS for several sites surrounding the vortex are 
shown in black, red, green, and blue. The orange line repre- 
sents TLDOS half-way from the vortex and the boundaries of 
the system. The black curve with a peak at i? = represents 
TLDOS at four sites surrounding the plaquette with the vor- 
tex. The parameters used in this figure are L — SOS, fi = 0.2t, 
VA ~ 0.25vf; the main structure of the graph is robust under 
the change of parameters of the model: in particular, at all 
cites near the vortex, except the four nearest neighbors, the 
TLDOS has additional peaks at energy Ac/3 — Ac/2. 



coherence peak energy. 



V. CONCLUSIONS 

We analyzed the properties of the mixed state in the 
tight-binding lattice c?-wave superconductors by consid- 
ering a quasiparticle spectrum of i) a vortex lattice and, 
separately, of ii) an isolated superconducting grain ac- 
commodating precisely a single vortex. To reduce the set 
of free parameters to a manageable number, we restricted 
ourselves to the simplest tight-binding model, described 
by only two parameters associated with the normal state 
(t, fj,), and additionally assume that the magnitude of the 
gap function A is spatially uniform. Within a set of sim- 
ple mean- field microscopic theories, such as those arising 
from an extended Hubbard ov at— J model, this assump- 
tion is well justified since for the parameters suited for 
to the cuprates, with their short coherence lengths, the 
self-consistent calculations show that the amplitude of 
the order parameter recovers its bulk value at distances 
of only a few lattice spacings. 

We find that the low-energy properties of the spectrum 
are described by Simon-Lee scaling only on average, and 
that both the energy dispersion and the (local) density of 
states experience oscillations as a function of the chemi- 
cal potential and a magnetic field. The magnitude of the 
oscillations in the energy levels behaves as as a func- 
tion of magnetic length, and therefore it is of the same 



order as the average Simon-Lee part of the dependence 
of Enk on I. We find that in all cases, these oscillations 
can be well described by a modified Simon-Lee scaling 
(|19|l . which includes additional 27r-periodic dependence 
of the energy levels on kijL. This modification is shown 
to be a consequence of the diverging solutions of the "lin- 
earized" Hamiltonian. A careful treatment of these diver- 
gencies is needed^i, lest one underestimates the quantita- 
tive importance of the inter-nodal and formally sublead- 
ing intra-nodal matrix elements of the Hamiltonian. As 
a result, the actual expansion parameters of the theory 
include kp^, rather than only kpl, and the scaling of the 
quasiparticle energy spectrum with respect to magnetic 
length I is consequently modified. 

In addition, we analyzed the large-energy, short- 
distance features of the quasiparticle spectrum. We 
found that - apart from the four sites surrounding a vor- 
tex where the thermally broadened TLDOS exhibits zero 
bias conductance peak (ZBCP), in agreement with Refs. 
00 - the TLDOS on all other sites in the vicinity of 
a vortex instead show peaks at sub-gap energies. The 
energy of these sub-gap peaks does not depend on mag- 
netic field, and is determined only by the parameters of 
the band structure and A. 

Finally, we examined how the TLDOS is modified 
when the amplitude of the d-wave gap function is var- 
ied locally in the vicinity of the vortex core and found 
that the suppression of the zero-bias peak corresponds to 
the enhancement of the d-wave gap function, which could 
arise through the locally enhanced^^'"^^ effective pairing 
interaction constant g. Such enhancement might result 
from the effect of the impurity atoms pinning the vortices 
or from a local distortion of atomic lattice by a vortex or 
fiuctuations of c?-wave order parameter. 
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APPENDIX: PHASE FACTORS IN THE 
TIGHT-BINDING HAMILTONIAN 

As mentioned in the main text, in principle the bond 
phase of the order parameter O^r' should be determined 
self-consistently. It is convenient, however, to adopt a 
synthetic approach and approximate 6rr' by using the 
known solution </)(r) of the continuum Ginzburg-Landau 
vortex-lattice problem, whose explicit form is given for 
example in Appendix A of Ref. |2ll One may set O^r' = 
0[(r + r')/2] or O^r' = (0(r) -H(/)(r'))/2. The latter, which 
will be predominantly used by us in this article, requires 
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FIG. 11; Left: Point group transformations: equivalent 
points in the unit cell, which are obtained by action of the 
group operation g from Table [TTI on the reference point 1, 
are shown as solid circles. Large open circles denote vor- 
tices. Right: The symmetry transformations of tight-binding 
Hamiltonian Htb- If tpr is an eigenstate of energy E, then 
the states ^'^ are also eigenfunctions of the Hamiltonian with 
the same energy, but different momenta in the Brillouin zone 
listed in the right column of the table. For brevity only a 
half of all transformations is listed. The remaining half is ob- 
tained by applying operator P^r = 7ri/'-r to each operation 
in the Table. Thus, overall there are 16 symmetry operations 
requiring that the spectrum is 16-fold symmetric as shown in 
the right panel of Fig. 



explanation since 9 must be defined modulo 27r, 
in the form above 9rr' is defined only modulo tt. 
accurate definition reads 



exp(i6'rr') — exp 



dr" 



while 
More 



(A.l) 



where the integral is over the bond connecting r and r'. 
It is easy to show that this definition is consistent in the 
sense that exp(i0rr') — cxp(z0r'r)i a-nd provided that (f){r) 



does not change by more than tt along a bond, the phase 
6'rr' defined in this way is indeed the "average" of (/)(r) 
and 0(r') in the sense that the "average" is understood 
as the closest to either 0(r) (or equivalently ?!>(r')) of the 
two possible choices. 

Note that unlike 6'rr' , phases 0^ and are not deter- 
mined by any self-consistent procedure, and merely serve 
as a technical device to recast the Hamiltonian in a pe- 
riodic form; we are therefore free to assign their values 
according to our convenience. Without lost of generality, 
we choose them by simply evaluating continuous func- 
tions (j)A (r) and 0b (r) from the previous section on sites 
of the tight-binding lattice. 

Given definition HA.1(I . coefficients V^,, V^, , and Arr' 
can be easily found from (0 and explicit expressions for 
4'A{B) from Appendix A of Ref. I21I 



_ pij^ Vc{r")-dr" 



A,B 



(A.2) 



62 . 



: (vA(r")-VB(r"))-cir" (A. 3) 

Since VA{r) and VB(r) are periodioi^iSi, clearly so are 

Vr"r' and Ar'- 

The phase factors possess a number of discrete sym- 
metries, which result in the symmetry of the dispersion 
shown in the right panel of Fig. ^ The full set of the sym- 
metry operations consists of operations g shown inlllland 
eight additional transformations Pg, where Ptpr = 7rV'-r 
is an inversion operator around a midpoint between two 
(arbitrary) vortices A and B. Each transformation in- 
volves a point group operation gr shown in Fig. 1111 
which may be followed by complex conjugation and mul- 
tiplication by Pauli matrix a^. Thus, overall there are 16 
distinct points in the Brillouin zone with identical set of 
energy eigenvalues. 

Additionally, if ipr is an eigenstate of energy E, then 
f2'01r an eigenstate of energy {—E) and the same mo- 
mentum k. Thus, at each point the spectrum is symmet- 
ric as a function of energy. 
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